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Frustrated quantum XXZ spin chains with the next-nearest-neighbor (NNN) couplings are typ- 
ically deterministic many-body systems exhibiting Gaussian orthogonal ensemble (GOE) spectral 
statistics. We investigate energy diffusion for these spin chains in the presence of a periodically 
oscillating magnetic field. Diffusion coefficients are found to obey the power law with respect to 
both the field strength and driving frequency with its power varying depending on the linear re- 
sponse and non-perturbative regimes. The widths of the linear response and the non-perturbative 
regimes depend on the strength of frustrations. We have also elucidated a mechanism for oscillation 
of energy diffusion in the case of weakened frustrations. 



I. INTRODUCTION 

There exists an accumulation of studies on quantum 
dynamics of classically chaotic systems , e.g., kicked rota- 
tors, kicked spin-tops, hydrogen atoms in time-dependent 
electric field, and the standard map model, to mention a 
fewi Quantum suppression of energy diffusion, dynam- 
ical localization and other signatures of quantum chaos 
are notable in these dynamics. However, most of the 
systems treated so far are confined to those with a few 
degrees-of-freedom, and little attention is paid to dy- 
namics of quantum many-body systemsSi^ whose adi- 
abatic energy levels are characterized by Gaussian or- 
thogonal ensemble (GOE) spectral statistics, i.e., by a 
hallmark of quantum chaos. While some important con- 
tributions^iSiL&SiiSiii are devoted to dynamics of a kind 
of many-body systems, those systems are actually de- 
scribed by the random-matrix models, and not by de- 
terministic quantum Hamiltonians. It is highly desirable 
to explore dynamical behaviors of deterministic quan- 
tum many-body systems exhibiting GOE or GUE spec- 
tral statistics. 

On the other hand, the frustrated quantum spin sys- 
tems have been receiving a wide attention, and we can 
find their realization in s — i antiferromagnetic chains 
Cu(ampy)Br2i£ and (N 2 H 5 )CuCl3r^ and in s = \ trian- 
gular antiferromagnetsiii The high-lying states of these 
quantum many-body systems deserve being studied in 
the context of "quantum chaos." The advantage of the 
frustrated quantum systems is that one can expect quan- 
tum chaotic behaviors appearing already in the low en- 
ergy region near the ground state^£*i£ From the view- 
point of real physics of condensed matters, novel features 
observed in the low-energy region are very important 
and welcome. Recalling that in most of deterministic 
Hamiltonian systems quantum chaotic behaviors appear 
in high-lying states, the role of frustration is essential in 
the study of quantum dynamics from the ground state of 



deterministic many-body systems with GOE or Gaussian 
unitary ensemble (GUE) level statistics. 

In this paper, we investigate dynamics of XXZ quan- 
tum spin chains which have antiferromagnetic exchange 
interactions for the nearest-neighbor (NN) and the next- 
nearest- neighbor (NNN) couplings. The NNN couplings 
cause the frustration, i.e., difficulty in achieving the 
ground state, thereby attributing a name of frustrated 
quantum spin chains to these systems. In fact, the level 
statistics of the NNN coupled XXZ spin chains without 
an applied magnetic field has been studied intensively 
in Refs. ITtIITsI and it has been shown that GOE behav- 
ior, which is typical of quantum chaos, appears already 
in the low energy region near the ground state^S* The 
ground-state phase diagram is shown in Ref. 00] for the 
NNN coupled XXZ spin chains without a magnetic field. 

A natural extension of the research is to investigate dy- 
namics of the frustrated quantum spin chains with an ap- 
plied periodically oscillating magnetic field. We calculate 
a time evolution of the system starting from their ground 
state and analyze the nature of energy diffusion. We shall 
numerically exhibit the time dependence of energy vari- 
ance, and show how the diffusion coefficients depend on 
the coupling constants, the anisotropy parameters, the 
magnetic field and the frequency of the field. Further- 
more, to compare with the energy diffusion in the case 
of weakened frustrations, we also investigate dynamics of 
the corresponding energy diffusion in XXZ spin chains 
with small NNN couplings. 

The organization of the paper is as follows: In Sec. [HJ 
we briefly describe a numerical approach to obtain the 
time evolution operator. In Sec. IIHI we shall show the 
time dependence of energy variance starting from the 
ground state of the many-body system and explain a way 
to evaluate diffusion coefficients. Section IIVI elucidates 
how diffusion coefficients depend on field strength and 
driving frequency. Here power laws are shown to exist 
in the linear response and non-perturbative regions. Sec- 
tion [V] is devoted to a mechanism of oscillation of energy 
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diffusion. Conclusions are given in Sec. |Vj 



II. NUMERICAL PROCEDURE 

We give Hamiltonian for the NN and NNN exchange- 
coupled spin chain on L sites with a time-periodic oscil- 
lating magnetic field as 



H(t)=H + H 1 (t), 



(1) 



where 



i=i 

L 

+ J 2 J2(S*S* +2 + S]S« +2 + AS?S] +2 ) 



3 = 1 



(2) 



(3) 



i=i 



Here, = {l/2)af and (crj, ctJ, <tJ) are the Pauli ma- 
trices on the jth site; the periodic boundary conditions 
(P. B. C.) are imposed. The magnetic field -BJ on jth 
site along the z axis is chosen to form a traveling wave: 



B](t) = B sin 



uit 



¥) 



(4) 



The period of Eq. (JH as well as Eq. (@J is T = 27r/w. 
Because of the coexisting spatial P. B. C, however, the 
effective period of the adiabatic energy spectra is given by 
T' = T / L = 2ir/(uL). In other words, the period of the 
Hamiltonian operator is T, and the spectral flow of the 
eigenvalues has the effective period T". This periodicity 
property comes from the traveling- wave form of Eq. 
and is advantageous for our getting a sufficient number 
of relevant data in each period T. 

When Ji > and J 2 > 0, the unperturbed Hamilto- 
nian Tig without coupling to the magnetic field is trans- 
lationally invariant and corresponds to a frustrated anti- 
ferromagnetic quantum spin model exhibiting GOE level 
statistics(i£*i& If J 2 = and Bq — 0, it describes an 
integrable and non-frustrated model. Before calculat- 
ing energy diffusion, we have to consider the symme- 
tries of the model. We divide the Hamiltonian matrix 
to some sectors which have the same quantum numbers. 
In the Hamiltonian Eq.Q), total S z (S* ot ) is conserved. 
The eigenstates with different Sf ot are uncorrelated. On 
the other hand, the non-uniform magnetic field breaks 
the translational symmetry, and leads to mixing between 
manifolds of different wave- number values. 



Before proceeding to consider the time evolution of a 
wave function, we should note: If we use the original 
Hamiltonian 7i(t) = Ho + Hi(t) as it stands, the mean 
level spacing of eigenvalues would change depending on 
J 2 , A, and Bq. To see a universal feature of the energy 
diffusion, it is essential to scale the Hamiltonian so that 
the full range of adiabatic energy eigenvalues becomes 
almost free from these parameters. Noting that this en- 
ergy range for the original Hamiltonian is of order of L 
when J i — J 2 — A — 1 , we define the scaled Hamiltonian 
H {t) ~ Hq + Hi (t) so that the full energy range equals L 
at t = 0, which will be used throughout in the text. The 
Schodinger equation is then given by 

ih^m)) =H{t)m)) = [#o+#i(*)]iv>(*)>- (5) 

The solution of Eq. © consists of a sequence of the 
infinitesimal processes as 



\i/j(t)) = U(t;t- At)U{t- At;t-2At) 
■■■ U(2At;At)U(At',0)\ip(0)). 



(6) 



The initial state |^>(0)) is taken to be the ground state, 
since our concern lies in the dynamical behaviors starting 
from the many-body ground state. To calculate a time 
evolution operator U(t + At; t) for each short time step 
At, we use the fourth-order decomposition formula for 
the exponential operator^ 



U(t + At;t) 



S(-ip 5 At/h, t 5 )S(-ip 4 At/h, ti) 
S(~i P2 At/h, t 2 )S(-i Pl At/h, ti), (7) 



where, 



l 'xH 1 (t)\ , TT . ( xHAt) 
b(x, t) = exp — - — exp(xii ) ex P 



(8) 



Here, tj's and pj's are the following: 

(pi +P2 + ■ ■ ■ + Vj-i + Pj/2)At, 



^ = t 



P = Pl=P2=Pi =P5, 

= 0.4144907717943757 • • • 
p 3 = 1 - 4p. 



(9) 



The numerical procedure based on the above decompo- 
sitions is quite effective when Hi (t) and Hq do not com- 
mute and each time step is very small. Our computation 
below is concerned mainly with the system of L = 10, 
whose S* ot = 1 manifold involves 210 levels. To check 
the validity of our assertion, some of the results will be 
compared to those for the system of L 
whose manifold involves 364 levels. 



14 and S t z ot = 4 



III. TIME DEPENDENCE OF ENERGY 
VARIANCE 

We calculate time evolution of the state and evaluate 
energy variances at each integer multiple of the effective 
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FIG. 1: (Color online) Time evolution of energy diffusion for 
(a) L — 10 and (b) L — 14. The parameters are the following: 
Ji = J 2 = 1.0, A = 0.3, Bo = 1.0. 



period T 1 = T/L = 2ir/(uiL). As mentioned already, 
we choose the ground state as an initial state, follow- 
ing the spirit of real physics of condensed matters. This 
viewpoint is in contrast to that of the random matrix 
models where initial states are chosen among high-lying 
oneS| 5 i 6 i 7 i 8 i 9 i io Consequently, the energy variance of our 
primary concern is the variance around the ground state 
energy Eq and is defined by 

5E(t) 2 = m)\[H(t)-E ] 2 \m)- (10) 

Time evolution of 8E(t) 2 is shown in Fig.^ The param- 
eters except for u> are fixed. The larger u is, the faster 
the energy diffusion grows, which is consistent with our 
expectations. The details will be explained in Sec. IIVI 
For wide parameter values of the next-nearest-neighbor 
(NNN) coupling J2 and exchange anisotropy A, the early 
stage of quantum dynamics becomes to show the normal 
diffusion in energy space, i.e., a linear growth of 5E(t) 2 
in time. While we proceed to investigate this normal 
diffusion process, energy variances will finally saturate 
because the system size we consider is finite. On the 
other hand, energy variances can also saturate because 
of another reason, i.e., the dynamical localization effect 
associated with a periodic perturbation. 

During the first period, SE(t) 2 shows a linear growth 
in time as shown in Fig. ^ ( a )- The range of the linear 
growth is not sufficiently wide because the number of 
levels is not large enough for L = 10. However, if the 
number of levels as well as the system size is increased, 
the length of a linear region may be elongated. In fact, 
the linear growth of 8E(t) 2 during the first period can be 
recognized more clearly for L = 14 than for L = 10 [see 
Fig-ffl(b)]. The diffusion coefficient has to be determined 
much earlier than the time where saturation begins. We 
determine the diffusion coefficient D from the fitting 



8E(t) 2 = Dt + const. 



(11) 




(b) J 2 =0.2, 

ri 



A' 
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to some data points around the largest slope in the first 
period, where the normal diffusion is expected. 
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FIG. 2: Driving frequency dependence of the diffusion coeffi- 
cients. The chained line and the solid line are just eye guides 
for D oc uj 13 with (5 = 1 and 2, respectively. The symbols 
(o) are the average of the diffusion coefficients calculated for 
several values of A (0.3 < A < 0.8). The parameters are the 
following: L = 10, Ji = 1.0; (a) J 2 = 1.0, (b) J 2 = 0.2. 



IV. DIFFUSION COEFFICIENTS: 
DEPENDENCE ON FIELD STRENGTH AND 
FREQUENCY 

Since the time evolution of our system starts from the 
ground state, we consider non-adiabatic regions where 
inter-level transitions frequently occur. In other words, 
we suppress a near-adiabatic or the so-called Landau- 
Zener (LZ) region where the driving frequency u> is much 
smaller than the mean level spacing divided by Planck 
constant. Because of a large energy gap between the 
ground and first excited states, the near-adiabatic region 
cannot result in the notable energy diffusion and will be 
left outside a scope of the present study. 

Beyond the LZ region, however, so long as the changing 
rate X of a perturbation parameter is not very large^ 
the diffusion coefficient can be calculated using the Kubo 
formula. We call such a parameter regime "linear re- 
sponse" regime. In the linear response regime, D oc X 2 
(See, e.g., Refs. and 0). When X is large, how- 
ever, the perturbation theory fails. We call such a pa- 
rameter regime "non-perturbative" regime. In the non- 
perturbative regime, the diffusion coefficient is smaller 
than that predicted by the Kubo formulaiS According 
to Ref. |3, D oc X 1 with 7 < 1 in the non-perturbative 
regime. We note that X oc Bqoj in this paper since the 
perturbation is given by Eq. Both Refs. and0 are 
based on the random matrix models, which are utterly 
different from our deterministic one. 

Numerical results of diffusion coefficients in Fig. [21 are 
almost consistent with the argument of Ref. [jj Diffusion 
coefficients as a function of u> are shown in Fig. [21 In 
Fig. |2Ia), where J2 = 1.0 (i.e., the fully-frustrated case), 
D is larger as Bq is larger for a fixed value of u>. In 
a small-w regime, D oc cu 13 with (3 = 2, though (3 > 2 
for small Bq. The latter is merely attributed to the fact 
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FIG. 3: Dependence of the diffusion coefficients on the product of field strength Bo and driving frequency u> for (a) L — 10 
and (b) L = 14. The symbols (o) are the average of the diffusion coefficient calculated for several values of A (0.3 < A < 0.8). 
The parameters are Ji = Ji = 1.0; for the inset, Ji = 1.0 and J2 = 0.2. The chained line and the solid line are just eye guides 
for D oc (Bqu)) 13 with (3 = 1 and 2, respectively. Some error bars are too short to see. 



that the perturbation is too small to observe a sufficient 
energy diffusion when both lo and Bq are small. In a 
large- lo regime, (3 = 1. Namely, we observe that (3 = 
2 in the linear response regime and (3 = 1 in the non- 
perturbative regime. In fact, for a large-w regime, the 
increase of energy variances per effective period hardly 
depend on ui by the time when 5E(t) 2 starts to decrease. 
This explains the observation that D oc lo@ with (3 = 1 in 
both Fig.[3fa) and Fig-E^b) . Let us represent the increase 
of energy variances per effective period as A(SE 2 ). From 
the definition of D, i.e. Eq. (JTTJ}, D oc A{6E 2 )/T'. If 
A(SE 2 ) is constant, D oc lo. 

On the other hand, in Fig. Efb) where J2 = 0.2 (i.e., 
a weakly- frustrated case), the region with (3 = 1 is ex- 
panding. For small Bq, (3 > 2 in a small-u regime is 
the same as in the case of J 2 = 1.0. For small Bo and 
around lo ~ 1, D seems to rather decrease than increase 
especially in the case of J2 = 0.2. Some kind of local- 
ization would have occurred in the very early stage of 
energy diffusion for large lo and small Bq, leading to the 
suppression of D. 

It is seen more clearly in Fig. [3] how the behavior of 
D changes between a linear response regime and a non- 
perturbative regime. The diffusion coefficient D obeys 
the power law D oc (Bqlo)P with its power (3 being two 
in the linear response regime and (3 = 1 in the non- 
perturbative regime. For small Bqlo, the power law seems 
to fail because of some finite-size effects. These universal 
feature is confirmed in systems of larger size. Actually, 
D obeys the power law better for L = 14 [Fig.^b)] than 
L = 10 [Fig.|3Ja)]. In addition, error bars are shorter for 



L = 14 than L = 10. Here, we have used the data of 
lo < 1. We cannot expect meaningful results in a large-w 
regime since, as mentioned above, energy diffusion is not 
normal there. 

Figure |31 suggests that the strength of frustration 
should affect the range of the linear response regime. The 
linear response regime is shorter for J2 = 0.2 than for 
J2 = 1.0, while the non-perturbative regime is larger for 
J2 = 0.2 than for J2 = 1.0. In fact, when J2 = (i.e. the 
integrable case), D oc (B^lo) 13 with (3 = 1 for almost all 
the data in the same range of Bqlo as that of Fig. |3| 



V. OSCILLATION OF ENERGY DIFFUSION IN 
WEAKLY-FRUSTRATED CASES 

We shall now proceed to investigate oscillations of dif- 
fusion which occur in the non-perturbative regime of a 
weakly-frustrated case. Figure Ela) shows an example 
of oscillatory diffusion for J 2 = 0.2, which is compared 
with a non-oscillatory diffusion for J2 = 1.0. The two 
examples have the same set of parameters except for J2. 
However, the cases of J 2 = 1-0 and J 2 = 0.2 are in the lin- 
ear response regime and in the non-perturbative regime, 
respectively. The variance for both cases shows normal 
diffusion at the very early stage of time evolution. For 
J 2 = 1-0, the energy variance seems to saturate after a 
normal diffusion time. On the contrary, the energy vari- 
ance for J2 = 0.2 shows large-amplitude oscillations. To 
investigate more details, we introduce another definition 
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FIG. 4: Examples for time evolution of energy variances : (a) 
SE(t) 2 and (b) SE(t) 2 (see text). Solid lines are for J 2 = 1.0; 
Broken lines, J 2 = 0.2. The parameters are the following: 
L = 10, Ji = 1.0, A = 0.3, B = 1.5, w = 0.5. 

of energy variance: 

sE(t) 2 = (m\[H(t) - (m\H{t)\m)] 2 \m)- (12) 

This follows a standard definition of the variance and 
quantifies the degree of energy diffusion around the time- 
dependent expectation of the energy Hamiltonian. The 
time evolutions of 8E(t) 2 corresponding to that of SE(t) 2 
are shown in Fig. |l£b). In the fully- frustrated case 
(J 2 = 1.0), the profile of 8E(t) 2 is similar to that of 
8E(t) 2 . This observation indicates that an occupation 
probability spread over the whole levels after normal dif- 
fusion of energy. 

On the contrary, in a weakly-frustrated case ( J2 = 0.2) 
in Fig. 21 8E(t) 2 shows small-amplitude oscillations re- 
flecting the large-amplitude oscillations of SE(t) 2 . Most 
part of 5E(t) 2 for J2 = 0.2 is smaller than that for 
J2 = 1.0. Furthermore, minima of SE(t) 2 come just be- 
fore minima and maxima of 5E(t) 2 . These observations 
indicates the following: an occupation probability, which 
is diffusing slowly, clustering around the expectation of 
energy oscillates together with the expectation in the en- 
ergy space. To make the picture of such behavior clearer, 
let us consider an occupation probability described by 

p t (E n ) = \(<p n \m)\ 2 , (13) 

where \<fr n ) is the nth excited eigenstate of Ho: 

H Q \(f> n } = E n \4> n ). (14) 

When t — , Pt(E n ) is given by the Kronecker delta: 
Po(E n ) = SE n ,E , where Eq is the energy of the ground 



FIG. 5: Parts of energy spectra depending on adiabatically 
fixed time t with < t < T/4. Effective period is luT' = 
27r/10. The parameters are the following: L = 10, J\ = 1.0, 
A = 0.3, Bo = 0.8; (a) J 2 = 1.0, (b) J 2 = 0.2. 



state. As t increases, Pt(E n ) forms a wave packet in en- 
ergy space and moves to higher levels. When the wave 
packet reaches some highest level, it reflects like a soli- 
ton and moves back to lower levels. Such behavior is 
repeated, although the wave packet of P t (E n ) broadens 
slowly. We have actually watched this soliton-like behav- 
ior of P t (E n ) in a form of an animation. 

The picture discussed above is also supported by the 
adiabatic energy spectra in Fig. EJ Figures a) andOJb) 
correspond to fully- and weakly-frustrated cases, respec- 
tively. Much more sharp avoided crossings appear in 
Fig- 13b) than Fig. 0a)- Some energy levels appear to be 
crossing, although they are very close and never crossing 
in fact. At a sharp-avoided-crossing point, Landau-Zener 
formula for two adjacent levels is applicable. Then the 
nonadiabatic transition leads to one-way transfer of a 
population from a level to its partner, failing to result in 
the energy diffusion. For small- J2, therefore, successive 
sharp avoided crossings can suppress diffusion of energy. 

We believe that large-amplitude oscillations of SE(t) 2 
should be one of characteristic features of the non- 
perturbative regime in this finite frustrated spin system. 
In fact, similar oscillations of energy variance are seen 
for large u> and large Bq even when Ji = 1.0 though the 
energy variance rapidly converges after one or two peri- 
ods. How long such oscillations continue should depend 
mainly on J 2. 

It is a notable fact that, common to both J2 = 1.0 
and J2 = 0.2, the level-spacing distributions in Fig. 
show GOE behavior. This GOE behavior in the adia- 
batic energy spectra appears for an arbitrary fixed time 
except for special points such as t = T = 2tt/lu. This fact 
suggests that dynamics can reveal some various generic 
features of quantum many-body systems which can never 
be explained by level statistics. The level-spacing distri- 
butions in Fig. convey another crucial fact: they have 
been calculated for low energy levels because our interest 
is in the low energy region around the ground state. We 
have confirmed that the level-spacing distributions for all 
energy levels in the inset is also described by GOE spec- 
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FIG. 6: (Color online) Level-spacing distributions at t = n/4 
for lowest 300 levels from the ground state (about 10% of 
all 3003 levels). Blue histogram is for J2 = 1.0; Red bars, 
J2 = 0.2; Solid curve, GOE spectral statistics. The other 
parameters are the following: L = 14, Sjf ot = 1, J\ = 1.0, A = 
0.3, B = 0.8. The inset is for all levels when J 2 = 1.0. The 
numerical methods to obt ain the level-spacing distributions 
are referred in Refs. 



VI. CONCLUSIONS 

We have explored the energy diffusion from the ground 
state in frustrated quantum XX Z spin chains under the 
applied oscillating magnetic field. In a wide parameter 
region of next-nearest-neighbor (NNN) coupling J2 and 
exchange anisotropy A, the diffusion is normal in the 
early stage of diffusion. Diffusion coefficients D obey 
the power law with respect to both the field strength 
and driving frequency with its power being two in the 
linear response regime and equal to unity in the non- 
perturbative regime. In the case of weakened frustra- 
tions with small- J2 we find oscillation of energy diffusion, 
which is attributed to a non-diffusive and ballistic nature 
of the underlying energy diffusion. In this way, the en- 
ergy diffusion reveals generic features of the frustrated 
quantum spin chains, which cannot be captured by the 
analysis of level statistics. 
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